Quantum phase transitions of the diluted 0(3) rotor model 
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We study the phase diagram and the quantum phase transitions of a site-diluted two-dimensional 
0(3) quantum rotor model by means of large-scale Monte-Carlo simulations. This system has two 
quantum phase transitions, a generic one for small dilutions, and a percolation transition across the 
lattice percolation threshold. We determine the critical behavior for both transitions and for the 
multicritical point that separates them. In contrast to the exotic scaling scenarios found in other 
random quantum systems, all these transitions are characterized by finite-disorder fixed points with 
power- law scaling. We relate our findings to a recent classification of phase transitions with quenched 
disorder according to the rare region dimensionality, and we discuss experiments in disordered 
quantum magnets. 
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I. INTRODUCTION 

Quantum phase transitions occur at zero tempera- 
ture when a nonthermal parameter like pressure, mag- 
netic field, or chemical composition is varied. In 
the presence of defects, impurities, and other kinds 
of quenched disorder, the interplay between dynamic 
quantum fluctuations and static disorder fluctuations 
can lead to a variety of unconventional phenomena. 
Experimental examples include quantum Ising spin 
glasses, 1,2 heavy-fermion intermetallic compounds 
and other itinerant quantum magnets^ as well as high- 
temperature superconductors)^ the metal-insulator 
transition in metal-oxide-semiconductor field effect tran- 
sistors (MOSFETs))i2*ii and superconductor- insulator 
transitions in thin films. 12 

Quenched disorder has interesting consequences al- 
ready at classical phase transitions. In the early years, 
it was thought that impurities always destroy a critical 
point, because the system divides itself up into spatial 
regions that undergo the transition at different temper- 
atures (see discussion in Ref. 13 and references therein). 
However, it soon became clear that in classical sys- 
tems with short-range disorder correlations, the transi- 
tion generically remains sharp. Harris^* derived a cri- 
terion for the stability of a clean critical point against 
disorder: If the correlation length exponent v fulfills the 
inequality v > 2/d, where d is the spatial dimensionality, 
the critical behavior is not influenced by weak disorder. 
Harris' idea can be generalized to form the basis of a 
classification of critical points according to the behav- 
ior of the disorder strength under coarse graining^ The 
first class contains systems that fulfill the Harris crite- 
rion. In these systems, the disorder decreases without 
limit under coarse graining. The critical behavior is gov- 
erned by a clean renormalization group fixed point, and 
macroscopic observables are self-averaging. The other 
two classes can occur when the clean system violates the 



Harris criterion. In systems belonging to the second class, 
the disorder strength approaches a nonzero constant for 
large length scales, corresponding to a fixed point with 
finite disorder. The critical exponent are thus different 
from that of the corresponding clean system, with the 
new dirty correlation length exponent fulfilling the in- 
equality v > 2/d. 16 Moreover, macroscopic observables 
are not self-averagingiiZii& Finally, the third class con- 
tains systems in which the disorder strength (counter- 
intuitively) increases without limit under coarse grain- 
ing. The resulting infinite-randomness fixed point has 
unconventional properties including exponential rather 
than power-law scaling and very broad distributions of 
macroscopic observables ,A2i22i 

At zero-temperature quantum phase transitions, order- 
parameter fluctuations in space and time must be 
considered. 21,22 Quenched disorder is time-independent, 
it is thus perfectly correlated in one of the relevant dimen- 
sions, the (imaginary) time dimension. Because these 
correlations increase the effects of the disorder, quan- 
tum phase transitions are generically more strongly af- 
fected by disorder than classical transitions, potentially 
resulting in unconventional behavior. One of the earliest 
explicit examples was the random transverse-field Ising 
chain 19 : 20 : 23 (or the equivalent McCoy- Wu model 2 ^* 2 ^). 
This system belongs to the third of the classes discussed 
above, i.e., the critical point is of infinite-randomness 
type. The dynamical scaling is activated with the corre- 
lation time £ T and correlation length £ being related by 
ln£ r ~ (In contrast, at conventional critical points, 
this relation is a power law, £ T ~ £ z , with a universal 
dynamical exponent z). Analogous behavior has been 
found, e.g., in the two-dimensional transverse- field Ising 
model 15,26 and in quantum Ising spin glasses i 27 : 28 

An important aspect of phase transitions in disordered 
systems are the so-called rare regions, large spatial re- 
gions that are devoid of impurities or more strongly cou- 
pled than the bulk system. These regions can be in the 
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ordered phase even though the bulk system is still in the 
disordered phase. Griffiths*^ showed that this leads to a 
singularity (the Griffiths singularity) in the free energy in 
an entire parameter region (the Griffiths region or Grif- 
fiths phased) close to the phase transition. In generic 
classical systems with short-range disorder correlations, 
thermodynamic Griffiths effects are very weak because 
the singularity in the free energy is only an essential one. 
They are therefore probably unobservable in experiment. 
However, disorder correlations can greatly enhance the 
rare region effects MtSk 

Since quenched disorder is perfectly correlated in the 
(imaginary) time direction, quantum phase transitions 
are expected to display stronger rare region effects than 
classical transitions. Indeed, in the above-mentioned 
random quantum Ising systems, the Griffiths singulari- 
ties are of power-law type with the susceptibility diverg- 
ing over a finite parameter ran g e ,i2iSLSiSL2£ In itinerant 
quantum magnets, rare region effects can be even more 
dramatic. For Ising symmetry, the sharp quantum phase 
transition is destroyed by smearing^ because sufficiently 
large rare regions stop tunneling. The same also hap- 
pens in classical Ising magnets with plane defects^2**S 
and at certain nonequilibrium phase transitions^**^ A 
recent review of these and other rare region effects can 
be found in Ref . 13(1 

In systems with continuous order parameter sym- 
metry, the situation is more complex. The ground 
states of certain one-dimensional quantum spin chains 
are controlled by infinite-randomness fixed points^ On 
the other hand, in dimensions d > 2, the stable low- 
energy fixed point of random Heisenberg models has been 
shown to be conventional^*^ Preliminary renormaliza- 
tion group results*^ for the critical point in these mod- 
els suggested that the infinite-randomness fixed point is 
unstable, implying more conventional behavior. This 
agrees with Monte-Carlo simulations of diluted single- 
layer^Mi or bilayer^S* 4 *^ quantum Heisenberg antiferro- 
magnets that did not show indications of exotic scaling. 
Note, however, that inhomogeneous bond disorder can 
induce a new quantum-disordered phase that can be un- 
derstood as a quantum Griffiths phase^ 4 * 4 *^ 

In this paper, we report the results of large-scale 
Monte-Carlo simulations of a diluted 0(3) quantum rotor 
model in two space dimensions. We find that the system 
has two quantum phase transitions, a generic one for di- 
lutions below the lattice percolation threshold p c and a 
percolation type transition right at p c . Both transitions 
and the multicritical point that separates them display 
conventional power-law critical behavior. For the generic 
transition, the critical exponents are universal, i.e., inde- 
pendent of the dilution. A short account of part of this 
work has already been published in Ref. |4y. The present 
paper is organized as follows: In section^ we introduce 
the model and summarize the scaling theories for conven- 
tional and infinite-randomness critical points. The simu- 
lation method and our results for the phase diagram and 
the critical behavior of the quantum phase transitions 
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FIG. 1: Phase diagram of the diluted bilayer Heisenberg 
antiferromagnet, as function of J±/J\\ and dilution p. The 
dashed line is the percolation threshold, the open dot is the 
multicritical point of Refs. E^E^ . The arrow indicates the 
QPT studied here. Inset: The model: Quantum spins (ar- 
rows) Sj,i and Si, 2 reside on the two parallel square lattices. 
The spins in each plane interact with the coupling strength 
J||. Interplane coupling is J±. Dilution is done by removing 
dimers. 

are presented in section lTTTl In the concluding section llVl 
we relate our results to a general classification 47 of dirty 
phase transitions, and we consider experiments. 

II. THEORY 

A. Diluted quantum rotor model 

We consider a site-diluted O(N) quantum rotor model 
defined on a square lattice. Its quantum Hamiltonian is 
given by2£ 

Hq = U e * U ~ J e * e J hi ' A i • w 

Here, is an ./V-component unit vector at site i. Con- 
jugate momenta pi are defined via the usual canoni- 
cal commutation relations [h a ,pP] — iS a p on each site 
i. (a,/3 = 1...N are the component indices, and we 
work in units in which Ti = 1.) The components of 
the angular momentum L of each rotor are given by 
L"/ 3 — fi a p 13 — hPp a . The site dilution is described by 
the independent random variable q which can take the 
value and 1 with probability p and 1 — p, respectively. 

Elementary quantum rotors do not exist in nature; 
rather, they arise as effective low-energy degrees of free- 
dom of correlated quantum systems. For example, 0(2) 
quantum rotor models describe superconducting Joseph- 
son junction arrays or bosons in optical lattices. An A^=3 
quantum rotor describes the states of an even number of 
antiferromagnetically coupled Heisenberg spins. A spe- 
cific example is provided by the bilayer quantum Heisen- 
berg antiferromagnet depicted in the inset of Fig. 2] This 
system is equivalent to an 0(3) quantum rotor model 
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FIG. 2: Sketch of the classical model @. The arrows rep- 
resent the classical spins and the tubes show the locations of 
the linear defects (vacancies). 

with each dimer (S^S^a) of spins at site i and layers 
1 and 2 being represented by a single rotor. The rotor 
coordinate corresponds to S^i — S^2 and the angular 
momentum L 7 ; corresponds to Sj.i + £5^2 (see, e.g., chap- 
ter 5 of Ref. [22). The 0(3) quantum rotors also describe 
double-layer quantum Hall ferromagnets. 

We now focus on the 0(3) site-diluted quantum rotor 
model on a square lattice. Since we will be mostly in- 
terested in the universal critical behavior, we map the 
quantum system onto a classical system in the same uni- 
versality class. This can be easily achieved via a path 
integral representation of the partition function^ The 
resulting classical system is a three-dimensional Heisen- 
berg model with the extra dimension representing the 
imaginary time coordinate of the quantum rotor model. 
Because the impurities in the quantum rotor model are 
quenched (i.e., time-independent), the defects in the clas- 
sical system are linear, i.e., the disorder is perfectly cor- 
related in the extra (imaginary time) direction (see Fig. 
|2J). Thus, our classical Hamiltonian reads: 

H = K ^2 e i e 3 n i,T ■ n i,r + K^2e t n iyT ■ n ijT+ i, (2) 

where rii. T is an 0(3) unit vector at the lattice site with 
spatial coordinate i and "imaginary time" coordinate r. 
The coupling constant f3K of the classical model is re- 
lated to the ratio U / J of the quantum rotor model. Here, 
j3 = 1/T where T is an effective "classical" temperature, 
not equal to the real temperature in the quantum model 
which is zero. We set K — 1 and drive the classical sys- 
tem through the transition by tuning the classical tem- 
perature T. 

As an aside, we note that in the above-mentioned bi- 
layer Heisenberg antiferromagnet the dilution has to be 
done by removing dimers of corresponding spins in the 



two layers because each dimer is described by a single 
rotor. In contrast, for site dilution, the physics changes 
completely: Random Berry phase terms with no classical 
analogue arise. They are equivalent to impurity-induced 
moments^ and those become weakly coupled via bulk 
excitations. Thus, for all dilutions below the percolation 
threshold, p < p c , the ground state shows long-range or- 
der, independent of the coupling constants! This effect is 
absent for dimer dilution, and both phases of the clean 
system survive for small dilution. 

The critical behavior of the Hamiltonian in the 
clean limit (p = 0) is in the three-dimensional classical 
Heisenberg universality class, with the correlation length 
exponent v w 0.7. It thus violates the Harris criterion 
v > 2/dj_ = 1 (dj_ = 2 because only dimensions in 
which there is randomness count for the Harris criterion) . 
Therefore, the critical behavior must change upon dilut- 
ing the lattice. 

B. Power-law vs. activated scaling 

In this subsection, we summarize the conventional 
and activated scaling scenarios at critical points with 
quenched disorder to the extent necessary for the analysis 
of our simulation results. 

At conventional (quantum) critical points, correlation 
length £ and correlation time £ T are related by a power- 
law, £ T oc £ z with z being the dynamical critical expo- 
nent. (Note that in the effective classical system ©, £ 
and £ T are the correlation lengths in the space-like and 
time-like directions, respectively.) This is referred to as 
power-law dynamical scaling. In contrast, at infinite- 
randomness critical points, the dynamical scaling is ac- 
tivated, i.e., the relation between correlation length and 
time is exponential, ln(£ T ) oc ^,19^ 

These differences in the dynamical scaling lead to anal- 
ogous differences in the finite-size scaling behavior of ob- 
servables. If we denote the linear system size in the two 
space-like dimensions by L and the size in the time-like 
dimension by L T , the finite-size scaling forms of the mag- 
netization per site to = |m| and the susceptibility % at a 
conventional critical point read 

to = L-^ l, fh c {tL 1 / u ,L T /L z ) , (3) 
X = L'<^ Xc (tL 1 ^,L T /L z ) . (4) 

Here, t is the dimensionless distance from the critical 
point; and /3, 7 and v are the critical exponents of mag- 
netization, susceptibility, and correlation length, respec- 
tively. At an infinite-randomness critical point, the scal- 
ing combination L T /L Z has to be replaced by ln(L T )/L^ 
leading to the finite-size scaling forms 

to = L- /3 / l/ TO A (i J L 1 ^,ln(L T )/^) . (5) 
X = LT</ v x A (tL 1 / v , ln(.L T )/ L^) . (6) 

In addition to magnetization and susceptibility we also 
calculate three quantities whose scale dimension is zero 
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which makes them particularly suitable for locating the 
critical point and extracting high precision values for the 
correlation length and dynamical exponents. The first 
such quantity is the Binder ratio. It is defined by 



where [••.]_„ denotes the disorder average and (...) de- 
notes the Monte-Carlo average for each sample. This 
quantity approaches well-known limits in both bulk 
phases (stable fixed points): In the ordered phase, all 
spins are correlated, and the magnetization has small 
fluctuations around a nonzero value. Therefore, (|m| 4 ) ~ 
(|m| 2 ) 2 , and the Binder ratio approaches 2/3. In the dis- 
ordered phase, the system consists of many independent 
fluctuators. Consequently, (|m| 4 ) can be decomposed 
using Wick's theorem. For 0(3) symmetry this gives 
(|m| 4 ) « (15/9)(|m| 2 ) 2 , and the Binder ratio approaches 
4/9. More generally, the Binder ratio is large if all spins 
are correlated and decreases if the system contains in- 
dependently fluctuation units. Because the Binder ratio 
has scale dimension zero, its finite-size scaling form is 
given by 

9av = ~g c (tl}l v ,L T lL z ) or (8) 
g av = ~g A {tL x ' v M{Lr)/L^) (9) 

for conventional scaling or for activated scaling, respec- 
tively. Two important characteristics follow from the 
scaling form and the discussion above* 4 ^^ (i) For fixed 
L, g av has a peak as a function of L T . The peak position 
L™ 1 marks the optimal sample shape, where the ratio 
L T /L roughly behaves like the corresponding ratio of the 
correlation lengths in time and space directions, £ T /£- (If 
the aspect ratio deviates from the optimal one, the sys- 
tem can be decomposed into independent units either in 
space or in time direction, and thus g av decreases.) At 
the critical temperature T c , the peak value <?^ ax is inde- 
pendent of L. Thus, for power law scaling, plotting g av 
vs. L T /L™ ax at T c should collapse the data, without the 
need for a value of z. In contrast, for activated scaling 
the g av data should collapse when plotted as a function of 
log(L r )/ log(L" lax ). (ii) For samples of the optimal shape 
(L T = L™ ax ), plots of g av vs. temperature for different 
L cross at T c . 

The other two quantities of scale dimension zero 
we consider, are the ratios of disconnected correlation 
lengths and system sizes in both space and time-like di- 
mensions. Here, the disconnected correlation lengths £ dis 
and arise from the disconnected correlation function 
(ni )T n_j- )T /). In contrast the usual, connected correlation 
lengths £ and £ r arise from the connected correlation 
function (n^n^,-') — (nj iT )(rij )T /). The finite-size scaling 
forms of our ratios for conventional and activated scaling 
read 

£ dis /L = X c {tL l l v ,L T /L z ) or (10) 
i dis /L = XAitL^MLr)/^) , (11) 



and 

£>/L T = Y c {tl}! v ,L T /L*) or (12) 
C S /L T = YaQL^MLt)/^) , (13) 

respectively. Calculating these quantities provides inde- 
pendent checks for the location of the critical point and 
for the exponents z (or ip) and v. 



III. SIMULATIONS 
A. Monte Carlo method 

In order to study the phase transitions of the effec- 
tive classical model (01 we perform large scale Monte- 
Carlo simulations. We use the efficient Wolff cluster 
algorithm 51 to reduce the effects of the critical slowing 
down close to the phase transition. This is possible be- 
cause the dilution-disorder does not introduce frustra- 
tion, and all interactions are ferromagnetic. We investi- 
gate linear sizes up to L = 120 in space direction and 
L T = 2560 in imaginary time direction, for impurity con- 
centrations p — i, i, |, | and p c — 0.407253 which is 
the lattice percolation threshold. For the larger dilutions, 
p = |, |, and p c we perform both Wolff and Metropolis 
sweeps to equilibrate small dangling clusters. 

The determination of averages, variances, and distri- 
bution functions of observables in disordered systems 
requires the simulation of many independent samples 
with different impurity configurations. Because of the 
huge computational effort involved^ one must carefully 
choose the number N$ of disorder realizations (i.e., sam- 
ples) and the number Ni of measurements during the 
simulation of each sample for optimal performance. As- 
suming full statistical independence between different 
measurements (quite possible with a cluster update), the 
variance o~\ of the final result (thermodynamically and 
disorder averaged) for a particular observable is given 

fc .,53.54 

of = (ol + <tVNi)/N s (14) 

where as is the disorder-induced variance between sam- 
ples and a i is the variance of measurements within each 
sample. Since the computational effort is roughly pro- 
portional to NiNs (neglecting equilibration for the mo- 
ment), it is then clear that the optimum value of Ni 
is very small. One might even be tempted to measure 
only once per sample. On the other hand, with too short 
measurement runs most computer time would be spent 
on equilibration. 

In order to balance these requirements we have used a 
large number N$ of disorder realizations, ranging from 
1000 to several 10000, depending on the system size and 
rather short runs of 100-200 Monte-Carlo sweeps, with 
measurements taken after every sweep. (A sweep is de- 
fined by a number of cluster flips so that the total number 
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FIG. 3: Phase diagram of the three-dimensional Heisenberg 
model J2J as function of temperature T and concentration p 
of linear defects. MCP is the multicritical point. The big dots 
mark the numerically determined transition points. The lines 
are guides for the eye. 



of flipped spins is equal to the number of sites, i.e., on 
the average each spin is flipped once per sweep.) The 
length of the equilibration period for each sample is also 
100 Monte-Carlo sweeps. The actual equilibration times 
have typically been of the order of 10-20 sweeps at maxi- 
mum. Thus, an equilibration period of 100 sweeps should 
be more than sufficient. 
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Binder ratio g a v as a function of L T 
Lower panel: Power-law scaling plot 
gavlg™v x vs. L T /L™ ax Inset: Activated scaling plot g a vfg™v X 
vs. y = log(L r )/log(L™ a:E ). The statistical errors of the data 
in the two main panels are smaller than the symbol size (see 
text). 



B. Results: phase diagram 

We start the discussion of our results by considering 
the phase diagram of the classical Hamiltonian in 
the dilution-temperature plane. To determine the criti- 
cal temperature T c for a given dilution p, we use a simple 
iterative procedure based on the properties of the Binder 
ratio g av discussed after (0. We start with a guess for 
the dynamical exponent z (or, alternatively ip for acti- 
vated scaling). We then perform a number of simulation 
runs to calculate g av as a function of temperature for 
samples whose linear sizes fulfill L T oc L z . The approxi- 
mate crossing of the g av vs. T curves for different L gives 
an estimate for T c . At this temperature, we now calcu- 
late g av as a function of L T for fixed L. The maxima of 
these curves give an improved estimate for the "optimal 
shapes", i.e., for the dynamical exponent z. This proce- 
dure can be repeated until the estimate for T c converges. 
Typically, only two to three iterations were necessary for 
the desired accuracy. 

The resulting phase diagram is shown in Fig. [3J As 
expected, T c (p) decreases with increasing dilution p, but 
the ordered phase survives up to the lattice percolation 
threshold p c with T* — T c (p c ) > 0. Thus, the classical 
Heisenberg model @ with linear defects has two phase 
transitions, viz., a generic transition for p < p c and a per- 
colation type transition for p = p c and T < T* . They are 
separated by a multicritical point at (p c ,T*). Analogous 
behavior was found in the dimer-diluted bilayer quantum 



Heisenberg antiferromagnet 4 -' 4 ^ (see also Fig.0. 

We have carried out detailed investigations of both 
transitions and of the multicritical point. Our results 
for the critical behaviors will be presented in the next 
three subsections. 



C. Generic critical point for p < p c 

To determine the critical behavior of the generic tran- 
sitions and to test its universality, we have considered 
four different impurity concentrations, p — i, g, |, and 
For each concentration we have performed two types 
of simulations: The first consists of runs right at T c (p) 
for systems of different sizes L and L T with varying as- 
pect ratio L T /L. The finite-size scaling properties of g av , 
m, and \ allow us to extract the dynamical exponent 
z, as well as (3/v and y/v. In the second set of sim- 
ulations, we vary the temperature over a range in the 
vicinity of T c , but we consider only samples of "optimal 
shape" , L T oc L z , using the value of z found in the first 
part. Finite-size scaling then yields the correlation length 
exponent v. 

Let us start by discussing the behavior of the Binder 
ratio g av right at the critical temperature T c (p). The up- 
per panel of Fig.^shows g av as a function of L T for vari- 
ous L = 5 to 100 and dilution p = \ at T = T c = 1.1955. 
The statistical error of g av is below 0.1% for the smaller 
sizes and not more than 0.2% for the largest systems. As 
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FIG. 5: L™ ax /L vs. L for dilutions p = ±, §, f and |. Solid 
lines: Fit to L™ a * = aL z (l + 6L"" 1 ) with z = 1.310(6) and 
lji = 0.48(3). The statistical errors of the data are well below 
a symbol size. 



FIG. 6: m vs. L at T c for optimally shaped samples at 
dilutions p = |, i, | and |. The statistical error of the 
data is well below a symbol size. Solid lines: Fit to m = 
aL~ p/v {l + bL-" 1 ) with 0/v = 0.53(3) and uji = 0.48. 



expected at T c , the maximum Binder ratio for each of 
the curves does not depend on L. We now discriminate 
between power-law dynamical scaling (JSJ) and activated 
dynamical scaling To this end, we plot g av /9al ax 

as a function of L T jL™ ax in the lower panel of Fig. 0] 
The data scale extremely well, giving statistical errors 
of i™ ax in the range between 0.3% and 1%. For com- 
parison, the inset shows a plot of g av as a function of 
log(L T )/ log(L™ ax ) corresponding to activated dynami- 
cal scaling ®. Plotted this way, the data clearly do not 
scale. The results for the other impurity concentrations 
p — |, | , | are completely analogous. 

This analysis establishes that the dynamical scaling at 
the generic transition is of conventional power-law type. 
We now proceed to determine the dynamical exponent z 
and to study whether or not it is universal, i.e., indepen- 
dent of the dilution p. According to (JSJ, the maximum 
positions L™ ax should depend on L via a power law with 
the exponent z. In Fig. we plot L™ ax vs. L for all 
four dilutions p. The curves show significant deviations 
from pure power-law behavior which can be attributed to 
corrections to scaling due to irrelevant operators. In such 
a situation, a direct power-law fit of the data will only 
yield effective exponents. To find the true asymptotic 
exponents we take the leading correction to scaling into 
account by using the ansatz L" lax (L) = aL z (l + bL^^ 1 ) 
with universal (dilution-independent) exponents z and 
lo\ but dilution-dependent a and b. A combined fit of all 
four curves gives z = 1.310(6) and u)\ — 0.48(3) where 
the number in brackets is the statistical error of the last 
given digit. The fit is of high quality (x 2 ~ 0.7) and 
robust against removing complete data sets or removing 
points form the lower or upper end of each set. We thus 
conclude that the asymptotic dynamical exponent z is 
indeed universal. Note that the leading corrections to 
scaling vanish very close to p — | ; the curvature of the 
L™ ax (I) curves in Fig. |S] is opposite above and below 
this concentration. A straight power-law fit of the L" lax 
vs. L curve for this dilution gives z = 1.303(3) in good 
agreement with the value from the global fit. 
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FIG. 7: x vs - L at T c for optimally shaped samples at 
dilutions p — |, |, | and ~. The statistical error of the 
data is well below a symbol size. Solid lines: Fit to x = 
aL 7/ "(l + bL'" 1 ) with j/v = 2.26(6) and wi = 0.48. 



We now turn to the exponents j3/v and 7/2A Accord- 
ing to eqs. © and 10}, they can be obtained from the 
L-dependence of the magnetization and susceptibility at 
T c of the optimally shaped samples (L T = L™ ax ). In 
Fig. we plot the magnetization data for all four di- 
lutions p — |, i, and i. The statistical error of m 
is below 0.5%. The plots do not show strong curvature, 
and straight power-law fits give values between 0.495 and 
0.568 for the exponent (3 /v. This could be taken as an in- 
dication for nonuniversal behavior. However, given that 
z is universal, we have also attempted a combined fit 
of all four curves to m(L) = cL _/3/,y (l + eLL~ w ) with 
universal (3/v and u> (fixed to the value lo\ = 0.48 65 ) 
but dilution-dependent c and d. The combined fit works 
well (x ~ 1-6) and gives an asymptotic exponent of 
f3/v = 0.53(3). For comparison, a straight power-law 
fit for p = j (which is the dilution where the correc- 
tions to scaling approximately vanish, see above) gives 
(3/v = 0.527 in good agreement with the value from the 
global fit. We thus conclude that the data display no 
indication of a nonuniversal (3 /v. 

The analogous plot for the susceptibility data is shown 
in Fig. [7| Here, the statistical error of the data is below 
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FIG. 8: Binder ratio g av (left) and £,f ts / ' L T (right) as functions 
of temperature T for dilution p = 0.2. System sizes range from 
L = 5 to L — 100 with increasing slope. 



1%. Straight power law fits give effective values between 
2.02 and 2.28 for j/u. The combined fit of all four curves 
to the ansatz x(L) = el?l v {\ + fL~ u ) with ui = 0.48 
gives a universal asymptotic exponent j/v = 2.26(6). 
Again, a straight power-law fit to the data for p = | 
gives a value (viz., j/v = 2.22) in good agreement with 
the global fit. 

The exponents (3/v, j/v and z are not all independent 
from each other; they must fulfill the hyperscaling rela- 
tion 2(3/1/ + 7/V = d + z. Our values (3/v = 0.53(3), 
7/V = 2.26(6), and z = 1.310(6) fulfill the hyperscaling 
relation within the error bars, indicating that they can 
indeed be asymptotic values rather than effective expo- 
nents. 

After having discussed the simulations right at T c , we 
now vary the temperature over a range in the vicinity 
of T Cl but we consider only samples of "optimal shape" 
L T = l™ ax to keep the second argument of the scaling 
functions in eqs. 10), I0J, ©, lfTtj|) . and (fTSf) constant. 
Fig. [S] shows the Binder ratio g av and the ratio / L T 
as functions of temperature for for various L = 5 to 100 
and dilution p = i. (The correlation lengths have been 
calculated in the usual way via the lowest Fourier com- 
ponents of the spin-spin correlation functionS^&) The 
Binder ratio shows a near perfect crossing point, i.e., the 
corrections to scaling for this quantity are very small. 
In contrast, / L T displays larger corrections to scaling 
indicated by the drift of the the crossing point between 
different / ' L T curves with L. The behavior of £ dls /L 
(not shown) is very similar to that of /L T . 

We have therefore used a scaling analysis of the Binder 
cumulant as our main tool for determining v. Fig. [5] 
shows a scaling plot of g av vs. (T — T c )xl for impurity 
concentration p = ^ . (Here xl is the scaling factor nec- 
essary to collapse the data onto a master curve.) The 
quality of the scaling is very good, comparable to that 
in Fig. 0] However, since the scaling function lacks the 
characteristic maximum, the error of the resulting scal- 
ing factor xl is somewhat larger (1 to 2%) than that of 
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FIG. 9: Scaling plot of g av vs. (T - T c )x L for p = 0.2. x L 
is the factor necessary to scale the data onto a master curve. 
The statistical errors of g av are well below a symbol size. 




FIG. 10: Scaling factor vs. L for four disorder concentrations 
p = |, |, | and |. The statistical errors of xl are smaller 

-l/uj 



then a symbol size. Solid lines: Fit to xl 
with v = 1.16(3) and lo 2 = 0.5(1). 



gL 1/l/ (l + hL- 



L™. The data for other dilutions p = g, | and | lead 
to analogous scaling plots. To determine the correlation 
length exponent ^, we plot the scaling factor Xl vs. L for 
all four dilutions in Fig. Similar to the L™ ax vs. L 
curves in Fig.0] the xl vs. L curves show significant cur- 
vatures necessitating an ansatz xl = cL 1 ' v {l + dL^^ 2 ) 
that includes corrections to scaling. A combined fit of all 
four curves to this ansatz (with universal v and LO2) gives 
v = 1.16(3) and u>2 = 0.5(1). As above, the fit is robust 
and of high quality (x 2 ~ 1.2). Importantly, as expected 
for the true asymptotic exponent, v fulfills the inequality 
v > 2/d±=lm& Note that both irrelevant exponents u>\ 
and Ct>2 agree within their error bars, suggesting that the 
same irrelevant operator controls the leading corrections 
to scaling for both z and v. For comparison, we have also 
performed a straight power-law fit for the dilution where 
the corrections to scaling approximately vanish, p = j. 
It gives v — 1.12(3) in agreement with the value from the 
global fit. 

We have also performed an analogous scaling analysis 
of S^ 18 /L T . Because of the larger corrections to scaling, 



8 



30 




FIG. 11: Power-law scaling plot g a v/gZ ax vs. L T /L™ ax for 
p — p c = 0.407253 and T = 0.5. Statistical errors range from 
3 * 10 -4 for the small sizes to about 6 * 10~ 4 for the larger 
systems, as indicated. 



FIG. 12: L™ ax /L vs. L for dilution p = p c = 0.407253 and 
T — 0.5. The statistical error of the data is about a symbol 
size. Solid line: Power-law fit giving z = 1.83(3). 



the errors of the scale factors (i.e., slopes of the curves) 
are significantly higher, but within the error bar, the 
value of v agrees with that determined from g av . 



D. Percolation transition at p = p c 

After having discussed the generic dirty quantum rotor 
phase transition realized for p < p c , we now turn to the 
percolation-type transition occurring for p — p c and T < 
T* . At this transition, the dynamical fluctuations of the 
rotors are noncritical; and the critical behavior is due to 
the critical geometry of the percolating lattice. Vojta and 
SchmalianSi have developed a complete scaling theory for 
this percolation quantum phase transition. They have 
also calculated exact exponent values for the case of two 
space dimensions, viz., z = 91/48, v = 4/3, /3 = 5/36, 
and 7 = 59/12. 

In this subsection, we test these theoretical predic- 
tions by performing simulations at p — p c — 0.407253 
and T — 0.5. For two reasons, these calculations re- 
quire significantly higher numerical effort than those in 
the last subsection: (i) Due to the large value of the 
dynamical exponent z, the "optimal" linear size L™ ax 
in the time-like direction increases very rapidly with L. 
For our largest L = 80, the optimal L T turns out to be 
jjnax _ 2030 leading to a system of 13 million spins, 
(ii) The very strong geometric fluctuations of the lattice 
at the percolation threshold lead to noisier data. Thus a 
larger number of disorder realizations has to be averaged. 

Figure ^] shows the resulting scaling plot for the 
Binder cumulant as a function of L T /L™ ax for systems of 
sizes L = 9 to 80. The data scale reasonably well, but the 
quality is clearly less than that of the corresponding plot 
for the generic transition (Fig. 0}. Moreover there seems 
to be a small systematic broadening of the domes with in- 
creasing L which is likely caused by finite-size corrections 
to the critical lattice percolation problem. The resulting 
values of L™ ax have statistical errors of about 5%. To 
determine the dynamical exponent z, we plot L™ ax vs. 



L in Fig. The data can be well fitted by a power 
law giving an exponent of z = 1.83(3). The remaining 
small difference to the theoretical value 91/48 « 1.89 
can probably be attributed to corrections scaling for our 
rather small L. Indeed, a fit of the h™ ax vs. L data 
to the ansatz L™ ax {L) = aL 91 / 48 (l + bL~") is almost 
indistinguishable from the power- law fit. 

In addition to the Binder ratio we have also analyzed 
magnetization m and susceptibility \ for the optimally 
shaped samples (L T — L" lax ). In analogy to Figs. [6] and 
the L-dependencies of m and \ give the exponents 
(3/v and 7/f, respectively. The m vs. L curve shows 
noticeable upward curvature; and while a power-law fit 
gives 0/u = 0.15(3), using the ansatz m(L) — aL 5 / 48 (l + 
dL~ w ) actually leads to a significantly better fit (lower 
X 2 ). For the susceptibility, a power-law fit of the x vs. L 
data gives 7/^ = 3.51(5) which is somewhat smaller than 
the theoretical value of 59/16 w 3.68. However, a fit to 
the ansatz x(L) = e£ 59 / 16 (l + fL~ u ) is of comparable 
quality. 



E. Multicritical point 

In this last subsection on results, we consider the mul- 
ticritical point (p c , T*) that separates the generic transi- 
tion from the percolation transition. In contrast to the 
percolation transition, no quantitative theoretical pre- 
dictions are available for the multicritical point. Wc 
have performed simulations at p = p c = 0.407253 and 
T = T* = 0.791 for L = 9 to 64. Because of the simulta- 
neous presence of critical geometric and dynamical fluc- 
tuations, the necessary number of disorder realizations is 
even larger than for the percolation transition. We have 
used between 10000 and 50000 realizations, depending on 
system size. 

Figure ^| shows the resulting scaling plot for the 
Binder cumulant as a function of L T jL™ ax . The data 
scale very well, giving statistical errors for L" lax of ap- 
proximately 2%. Note that the data show a very slight 
systematic broadening of the domes with increasing L. 
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FIG. 13: Power-law scaling plot g a v/gTv X vs. L T /L™ ax for 
the multicritical point at p = p c = 0.407253 and T — T* — 
0.791. The statistical error of the data is about a symbol size. 



Exponent Generic Multicritical Percolation Perc. (theory) 



z 


1.310(6) 


1.54(2) 


1.83(3) 


91/48 


P/v 


0.53(3) 


0.40(3) 


0.15(3) 


5/48 


l/v 


2.26(6) 


2.71(3) 


3.51(5) 


59/16 


V 


1.16(3) 









TABLE I: Numerical results for the critical exponents. For 
the generic transition, the values are from the combined fit of 
all four datasets for dilutions 1/8, 1/5, 2/7, and 1/3. For the 
multicritical point and the percolation transition the values 
are from straight power-law fits. The numbers in brackets give 
the statistical error of the last given digits. For the percolation 
transition we also give the theoretical results from Ref. I5H . 
The numerical data are compatible with these values if one 
allows for corrections to scaling. 




FIG. 14: L™ ax /L vs. L for the multicritical point p = p c — 
0.407253 and T = T* = 0.791. The statistical error is well 
below a symbol size. Solid line: Power-law fit giving 2 = 
1.54(2). 



It is much weaker than for the percolation transition at 
T = 0.50 (Fig.HTjl. but probably, it can also be attributed 
to finite-size corrections to the critical lattice percolation 
problem. Fig. ^] shows a log-log plot of L™ ax vs L. 
The curve does not show a discernable deviation from 
a straight line, and a power-law fit gives z — 1.54(2). 
Given the fact that the same analysis gave a very slightly 
too small z-value at the percolation transition, the true 
asymptotic exponent may be a few percent higher than 
the fit result. We have also determined the exponents 
(3/v and "f/v from the L-dependencies of the magnetiza- 
tion and susceptibility for the optimally shaped samples. 
As at the percolation transition, the m(L) curve shows 
some upward curvature; and while a power-law fit gives 
P/u = 0.40(3), the true asymptotic exponent may be a 
bit lower. Unfortunately, our L-range is not wide enough 
for a stable fit to an ansatz that includes corrections to 
scaling (and thus has two unknown exponents and two 
prefactors). In contrast, the curve does not show 

any deviations from power-law behavior, and a fit gives 
7/1/ = 2.71(3). Again, from the analogy with the perco- 
lation transition, the true asymptotic exponent may be 
slightly higher. 



IV. CONCLUSIONS 

We have investigated the quantum phase transitions 
of a two-dimensional site-diluted 0(3) rotor model by 
performing large-scale Monte-Carlo simulations of the 
equivalent classical model, a three-dimensional classical 
Heisenberg model with linear defects. In this final sec- 
tion we summarize the results and relate them to a recent 
classification^^! of phase transitions with quenched dis- 
order. We also compare our findings with previous work 
on this and related problems; and we consider experi- 
ments. 

The two-dimensional site-diluted 0(3) rotor model has 
two quantum phase transitions, (i) a generic transition 
for dilutions p below the percolation threshold p c of the 
lattice and (ii) a quantum percolation transition at p c . 
These transitions are separated by a multicritical point. 
Our calculations have shown that the critical behavior of 
all these transitions is of conventional power-law type. In 
contrast, the Ising version of our model, the diluted 2d 
random transverse-field Ising model, shows an infinite- 
randomness critical point li^SS To study the generic 
transition, we have considered four different dilutions. 
Combined fits of all data sets have allowed us to system- 
atically include corrections to scaling. In this way, we 
have provided strong evidence that the critical behavior 
of the generic transition is universal, i.e., independent 
of the dilution. Because of the high numerical effort, 
we have considered only one data set for the percolation 
transition. This one data set does not allow us to extract 
the leading exponents and the corrections to scaling from 
a free fit of the numerical data. However, including cor- 
rections to scaling when fitting the data to the theory of 
Ref. |H3 leads to a very good agreement. For the multicrit- 
ical point, there is only one data set, and no quantitative 
theoretical results exist. Therefore, our results for the 
multicritical behavior are on somewhat less firm ground 
because they do not contain corrections to scaling. Our 
exponent values are summarized in table [I] 

Recently, a general classification has been suggested for 
phase transitions with weak (random- T c type) quenched 
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disorder and short-range interactionsi^iiH According to 
this classification, the type of critical behavior depends 
on the effective dimensionality cIrr of the defects or, 
equivalently, the rare regions. Three cases can be distin- 
guished. (A) If cIrr is below the lower critical dimension 
d~ of the problem, rare region effects are exponentially 
small. As a result, the transition is sharp, and the crit- 
ical point is of conventional power-law type. (B) In the 
second class, with dun = d~, rare regions are much more 
important. A sharp transition still exists, but the crit- 
ical point is controlled by an infinite-randomness fixed 
point with activated scaling. In addition, there are strong 
power-law Griffiths effects. (C) Finally, for djm > d~, 
the rare regions can order independently leading to a de- 
struction of the sharp phase transition by smearing. 

In the problem considered here, dun = 1 because 
the defects are linear. The lower critical dimension of 
the Heisenberg universality class is d~ — 2. Therefore, 
dun < d~, and our model should be in class A with 
conventional power-law critical behavior. Our numerical 
results are thus in complete agreement with the above 
general rare-region based classification scheme. 

Let us compare our results to previous work. The 
qualitative structure of the phase diagram, viz., the fact 
that long-range order survives for all dilutions up to and 
including the percolation threshold agrees with earlier 
quantum Monte-Carlo simulations for the bilayer quan- 
tum Heisenberg antifcrromagne^Si^ and with analytical 
results for diluted magnets 5 ^ as well as 0(2) rotors^ 
Sandvik^ and Vajk and Greveri 4 ^ studied the multicriti- 
cal point at p — p c . They found a dynamical exponent of 
z ?s 1.3 significantly lower than our result of 1.54. More 
recently, SandvikSi reported a somewhat larger value of 
z = 1.36, but it is still well below our result. The reasons 
for this discrepancy are presently not fully understood. 
Possible explanations include a failure of the quantum to 
classical mapping (which we consider unlikely) and cor- 
rections to scaling of the lattice percolation problem. In 
this context it is worth noting that the scaling properties 
of the lattice enter our calculations in a different way 
than that of Ref. 0. Our analysis of the Binder ratio 
works with linear extensions in space an time directions, 
directly giving z. In contrast, Ref. l6ll analyzes the tem- 
perature dependence of the susceptibility and effectively 
measures Df/z with Df being the fractal dimension of 
the percolation cluster. It is clear that corrections to 
scaling, if any, will enter the two calculations very differ- 
ently. 

Vajk and Greveni^ also quoted exponents for p < p c . 
At dilution p = 0.25 they find z = 1.07 and v = 0.89, 
different from our results. However, as the authors of Ref. 
l43l pointed out, a value of v < 1 violates the inequality 
v > 2/d, indicating that it represents an effective rather 
than an asymptotic exponent. It would also be useful to 
compare our exponents with analytical results. To the 
best of our knowledge, the only quantitative result for 
the generic transition is a resummation of the 2-loop e- 
expansioniSS The predicted exponents significantly differ 



from ours; but they also violate the inequality v > 2/d, 
casting doubt on their validity. 

While no quantitative analytical results exist for the 
multicritical point, there is a complete scaling theory for 
the percolation transition in the diluted rotor model, 57 
and the exponents in two dimensions are known exactly. 
Our numerical data are in excellent agreement with the 
exact values if one allows for corrections to scaling. Even 
if corrections to scaling are not included, the differences 
between the theoretical and numerical exponent values 
are only a few percent. Originally, this scaling theory 
was thought to apply not only to the rotor model but also 
to a site-diluted Heisenberg antiferromagnet; and there 
is some numerical evidence 44 in support of the relation 
z = Df predicted by the scaling theory. However, a 
recent exact diagonalization and quantum Monte-Carlo 
studj«Si finds a dynamical exponent z w 2Df for the 
site-diluted Heisenberg antiferromagnet but z ~ Df for 
the dimer-diluted bilayer (which is equivalent to a rotor 
model) . 

Finally, we discuss experiments. Chemical doping, i.e., 
random replacement of magnetic by non-magnetic ions, 
e.g., Cu by Zn in YBa2Cu30g, in both single-layer and 
bilayer antiferromagnets realizes site rather than dimer 
dilution. As discussed at the end of section III Al this 
leads to random Berry phases and a completely differ- 
ent physical picture. The most promising way to achieve 
dimer dilution is the introduction of strong antifcrromag- 
netic intra-dimer bonds at random locations. Thus we 
propose to study magnetic transitions in bond-disordered 
systems; those transitions can be expected to be in the 
same universality class as the one studied here. One can- 
didate material - albeit 3d - is (Tl,K)CuCi3^ under pres- 
sure; interesting quasi-2d compounds are SrCu2 (603)2 
or BaCuSi2 0g, where suitable dopants remain to be 
found. Our results may also be interesting for some 
single-layer Zn doped cuprate antiferromagnets that have 
been speculated to be parametrically close to the multi- 
critical point of the rotor model42^ Moreover, the quali- 
tative properties of the phase diagram and the critical be- 
havior are also important for disordered Josephson junc- 
tion arrays or diluted bosons in optical lattices. 
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